#define PERFORMANCE_TESTING  0

//
//  Standard headers
// 
#include <iterator>   // for stream_iterators, inserter
#include <iostream>
	using std::cout;
	using std::cerr;
	using std::cin;
	using std::ofstream;
	using std::ifstream;
#include <unordered_map>
	using std::unordered_map;
#include <unordered_set>
	using std::unordered_set;
#include <set>
	using std::set;
#include <string>
	using std::string;
#include <vector>
	using std::vector;
#include <map>
	using std::map;
#include <utility>
#include <algorithm>
	using std::min;
#include <numeric>
#include <stdexcept>
	using std::runtime_error;
#include <limits>
	using std::numeric_limits;
#include <algorithm>
	using std::max;
#include <tuple>
	using std::tuple;
	using std::make_tuple;
	using std::get;
	using std::tie;

using std::make_pair;

//
//  My common headers
// 
#include <progress_indicator.h>
#include <commify.h>

//
//  Boost headers
// 
#include <boost/range/adaptor/map.hpp>
	using namespace boost::adaptors;
#include <boost/range/algorithm/max_element.hpp>
	using boost::max_element;
#include <boost/range/algorithm_ext/push_back.hpp>
    using boost::push_back;
#include <boost/algorithm/string/join.hpp>
#include <boost/algorithm/string/replace.hpp>
	using boost::replace_all;
#include <boost/shared_ptr.hpp>
	using boost::shared_ptr;
#include <boost/foreach.hpp>
#include <boost/range/algorithm.hpp>
	using boost::sort;
#include <boost/format.hpp>
		using boost::format;
		using boost::io::group;
#include <boost/algorithm/string/replace.hpp>
		using boost::algorithm::replace_all;
#include <boost/iostreams/device/file.hpp>
#include <boost/iostreams/stream.hpp>
#include <boost/iostreams/device/mapped_file.hpp>
		namespace io = boost::iostreams;

#include <boost/assign/list_of.hpp> // for 'list_of()'
		using namespace boost::assign; // bring 'list_of()' into scope

#include <boost/assert.hpp> 
#include <boost/lexical_cast.hpp> 
#include <boost/spirit/include/karma.hpp>
#include <boost/regex.hpp>

#include <boost/random/mersenne_twister.hpp>
	using boost::mt19937;
#include <boost/random/uniform_int.hpp>
	using boost::uniform_int;
#include <boost/random/variate_generator.hpp>	
	using boost::variate_generator;




	
		
#include "tbb/task_scheduler_init.h"
    // tbb::task_scheduler_init

#include "intervals_per_contig.h"
	using intervals_per_contig::t_interval;
	using intervals_per_contig::t_intervals;
	using intervals_per_contig::t_intervals_per_contig;
	using intervals_per_contig::t_vec_intervals;
	using intervals_per_contig::t_map_str_intervals;
#include "read_options.h"
#include "parse_contig_pos.h"

#include "intervals_per_contig.h"
	using intervals_per_contig::t_interval;
	using intervals_per_contig::t_intervals;
	using intervals_per_contig::t_intervals_per_contig;
	using intervals_per_contig::t_vec_intervals;
	using intervals_per_contig::t_map_str_intervals;
	using intervals_per_contig::t_interval_search_cache;


typedef unordered_map<string, long> 		t_contig_genome_offsets;










//!_______________________________________________________________________________________
//!     
//!		get_intervals_from_all_files
//!			converts from intervals per contig to genomic coordinates
//!     
//!		\author lg (1/11/2011)
//!     
//!     \param opt 
//!     \param intervals 
//!     \param contig_genome_offsets 
//!_______________________________________________________________________________________
void get_intervals_from_all_files(prog_options& opt, vector<t_intervals>& intervals, t_contig_genome_offsets& contig_genome_offsets)
{
	for (unsigned i = 0; i < opt.interval_file_names.size(); ++i)
	{

		//
		//	Read intervals for each file
		// 
		string file_name = opt.interval_file_names[i];
		opt.both() << "Reading intervals from " << file_name << "...\n";
		io::stream<io::mapped_file_source> f_intervals(file_name);
		if (!f_intervals)
			throw std::logic_error("Can't open " + file_name);
		t_intervals_per_contig temp_intervals_per_contig;
		read_intervals(	f_intervals, temp_intervals_per_contig, opt.zero_based_intervals /*one_based_coordinates?*/);

		// 
		// 	Convert	from per contig coordinates to per genome coordinates
		intervals.push_back(t_intervals());
		t_vec_intervals& curr = intervals.back().data;
		BOOST_FOREACH(t_map_str_intervals::value_type& k_v, temp_intervals_per_contig.data) 
		{
			long long contig_pos_to_genome = contig_genome_offsets[k_v.first];
			BOOST_FOREACH(t_interval& interval, k_v.second.data) 
				curr.push_back(t_interval(	interval.beg + contig_pos_to_genome,
											interval.end + contig_pos_to_genome));
		}


		opt.both() << "    " << commify(curr.size()) << " intervals\n";
		intervals.back().combine_overlapping();
		opt.both() << "    " << commify(curr.size()) << " non-overlapping intervals\n";
		opt.both() << "    " << commify(intervals.back().len()) << " bp\n";
	}
	opt.verbose(1) << "Finished reading intervals...\n";

}

//!_______________________________________________________________________________________
//!     
//!		get_offsets_into_genome_pos_per_contig
//! 		Saves the offsets required per contig to convert contig coordinates to
//! 		genomic ones.
//!     
//!		\author lg (1/11/2011)
//!     
//!     \param opt 
//!     \param contig_genome_offsets 
//!     
//!     \return int total genome size
//!_______________________________________________________________________________________

struct is_canonical_contig
{
	bool operator()(const string& contig)
	{
		static boost::regex canonical_only("([0-9]+)|X|x");
		return regex_match(contig, canonical_only);
	}
};
long long get_offsets_into_genome_pos_per_contig(prog_options& opt, t_contig_genome_offsets& contig_genome_offsets)
{
	unordered_map<string, t_interval> contig_sizes;

	string file_name = opt.contig_sizes_file_name;
	opt.both() << "Reading contig sizes  from " << file_name << "...\n";

	io::stream<io::mapped_file_source> f_intervals(file_name);
	if (!f_intervals)
		throw std::logic_error("Can't open " + file_name);
	t_intervals_per_contig sizes_per_contig;
	read_intervals(	f_intervals, sizes_per_contig, opt.zero_based_intervals /*one_based_coordinates?*/);
	vector<string> contig_names = sizes_per_contig.get_sorted_contigs();


	unsigned cnt_sizes = sizes_per_contig.size();
	if (cnt_sizes != contig_names.size())
		opt.both() << "WARNING! Some contigs have more than one size\n";

	//	only keep 1-99,X
	//contig_names.erase(std::partition(
	//								contig_names.begin(), contig_names.end(), is_canonical_contig()), 
	//					contig_names.end());
	opt.both() << "    " << contig_names.size() << " contigs\n";
	opt.log() << "    contigs = [ " << boost::join(contig_names, ", ") << " ]\n";


	// only take the first size
	long long genomic_offset = 0;
	BOOST_FOREACH(const string& contig_name, contig_names) 
	{
		t_interval interval = sizes_per_contig.data[contig_name].data[0];
		genomic_offset -= interval.beg;
		contig_genome_offsets[contig_name] = genomic_offset;
		//opt.verbose() << boost::format("contig =%1$14s, beg = %2$15s, end  = %3$15s, offset = %4$15s\n") 
		//								% contig_name
		//								% commify(interval.beg)
		//								% commify(interval.end)
		//								% commify(genomic_offset);
		genomic_offset += interval.end;
	}
	return genomic_offset;
}







//!_______________________________________________________________________________________
//!     
//!		get_category_name_index
//! 		returns the index of the category name and whether this category is newly
//! 		created
//!     
//!		\author lg (3/15/2011)
//!     
//!     \param category_name 
//!     \param category_names_to_indices 
//!     
//!     \return tuple<unsigned,bool> 
//!_______________________________________________________________________________________


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// 
//      Main
// 
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888 
int main (int argc, char *argv[])
{
	try
	{
		// 
		//  Read options
		//  	return if user chose help or version
		// 
		prog_options opt;
		if (!read_options(opt, argc, argv))
			return 0;




		//
		//	Get contig sizes 
		// 
		t_contig_genome_offsets contig_genome_offsets;
		unsigned genome_size = get_offsets_into_genome_pos_per_contig(opt, contig_genome_offsets);



		//
		//  Read intervals from all files
		// 
		vector<t_intervals> intervals;
		get_intervals_from_all_files(opt, intervals, contig_genome_offsets);

		// 
		// add whole genome to intervals as catchall
		//
		intervals.push_back(t_intervals());
		intervals.back().data.push_back(t_interval(0, numeric_limits<unsigned>::max()));
		opt.interval_type_names.push_back(opt.default_interval_type_name);

		// index interval type names
		typedef unordered_map<string, int> t_map_str_to_index; 
		t_map_str_to_index	interval_type_names_to_indices;
		BOOST_FOREACH(const string& n, opt.interval_type_names) 
			if (!interval_type_names_to_indices.count(n))
			{
				auto index = interval_type_names_to_indices.size();
				interval_type_names_to_indices[n] = index;
			}



		// 
		//	input files
		// 
		ifstream input_stream;
		//io::stream<io::mapped_file_source> input_stream;
		input_stream.open(opt.input_file_name.c_str());
		if (!input_stream)
			throw std::logic_error("Can't open " + opt.input_file_name);
		opt.both() << "Input= " << opt.input_file_name << std::endl;

			
		// 
		// Start parsing
		// 
		opt.verbose() << "Parsing lines..." << std::endl;
		string line;
		line.reserve(5000);

		unsigned curr_line_number = 0;
		unsigned cnt_sites = 0;
		auto grammar = create_contig_pos_grammar(curr_line_number);
		// progress
		t_progress_indicator progress(opt.verbose(2), 200000U, 29000000U);
		progress.use_timer();
		t_contig_pos res;
		
		// maximum offset for any contig
		long max_contig_genome_offsets =  *max_element(contig_genome_offsets | map_values);
		// using this offset for any unknown contig
		const long max_contig_length = 1000000000;


		// search cache for each interval type
		vector<t_interval_search_cache> search_cache(intervals.size());


		//
		//  interval type name indices to save
		// 
		set<unsigned> to_save_interval_type_indices;
		BOOST_FOREACH(const string& n, opt.to_save_interval_type_names) 
			to_save_interval_type_indices.insert(interval_type_names_to_indices[n]);

		// 
		// category names to save
		// 
		set<string> to_save_categories(opt.to_save_categories.begin(), opt.to_save_categories.end());

		//
		//	interval type indices / category to save
		// 
		set<tuple<unsigned, string> > to_save_interval_type_category_pairs;
		BOOST_FOREACH(auto interval_cat, opt.to_save_interval_type_name_category_pairs) 
			to_save_interval_type_category_pairs.insert(make_tuple(interval_type_names_to_indices[get<0>(interval_cat)], get<1>(interval_cat)));
		

		vector<shared_ptr<ofstream>>						output_streams;
		map<string, unsigned>   							output_file_name_to_stream_index;
		map<tuple<unsigned, string>, shared_ptr<ofstream>>	interval_category_to_output_stream;
		


		while (getline(input_stream, line))
		{
			++curr_line_number;

			// ignore empty lines or comments
			if (!line.length() || line[0] == '#')
				continue;

			res.chromosome.clear(); 
			res.category.clear();
			auto beg = line.cbegin();
			if (!parse_contig_pos(grammar, beg, line.end(), res))
				throw std::runtime_error((boost::format("Failed on line#%1$") % curr_line_number).str());
	


			// 
			// Convert to genome wide coordinates
			// 	handling unknown contigs
			// 
			auto ii = contig_genome_offsets.find(res.chromosome);
			if (ii == contig_genome_offsets.end())
			{
				max_contig_genome_offsets += max_contig_length;
				contig_genome_offsets[res.chromosome] = max_contig_genome_offsets;
				ii = contig_genome_offsets.find(res.chromosome);
			}
			long genome_pos = res.pos + ii->second;
	
			//
			//	get interval type 
			// 
			unsigned interval_type = intervals.size() - 1;
			for (unsigned test_interval_type = 0; test_interval_type < intervals.size(); ++test_interval_type)
				if (intervals[test_interval_type].contains(genome_pos, search_cache[test_interval_type]))
				{
					interval_type = test_interval_type;
					break;
				}

			auto interval_category = make_tuple(interval_type, res.category);

			//
			//	Not yet saved to any stream
			// 
			if (!interval_category_to_output_stream.count(interval_category))
			{
				// 
				//	If don't need to save, continue
				// 
				if (!to_save_interval_type_category_pairs.count(interval_category)	&&
					!to_save_categories.count(res.category)							&&
					!to_save_interval_type_indices.count(interval_type))
					continue;

				//
				//	get stream name
				// 
				string output_file_name = opt.output_file_name;
				replace_all(output_file_name, "[INTERVAL_TYPE]", opt.interval_type_names[interval_type]);
				replace_all(output_file_name, "[CATEGORY]", res.category);

				//
				//	output stream is already there, use that
				// 
				if (output_file_name_to_stream_index.count(output_file_name))
				{
					unsigned strm_index = output_file_name_to_stream_index[output_file_name];
					interval_category_to_output_stream[interval_category] = output_streams[strm_index];
				}
				else
				{
					// 	create new stream
					//  link file name to that
					//  link interval / category to stream
					output_streams.emplace_back(new ofstream(output_file_name.c_str()));
					output_file_name_to_stream_index[output_file_name] = output_streams.size();
					interval_category_to_output_stream[interval_category] = output_streams.back();

				}
			}


			if (opt.add_interval_to_output)
				*interval_category_to_output_stream[interval_category] 
									<< line << "\t" << opt.interval_type_names[interval_type] << "\n";
			else
				*interval_category_to_output_stream[interval_category] << line << "\n";


	
			cnt_sites += 1;

		}



	}
	catch(const std::exception& e)
	{
		// indent error string
		string error_string = e.what();
		replace_all(error_string, "\n", "\n\t");
		std::cerr
				<< "\n"<< string(60, '=') << "\n"
				<< "\nError:\n\t" << error_string 
				<< "\n\n" << string(60, '=') << "\n\n";
		return -1;
	}
	return 0;
}


